library(ggplot2)
  library(reshape2)
  library(dplyr)
  library(tidyr)
  library(GGally)
  library(grid)
  "%&%" = function(a,b) paste(a,b,sep="")
  source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan-Paper/scripts/Heather/make-figures/multiplot.R')
  my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
  fig.dir <- '~/GitHub/GenArch/GenArchPaper/Figures/'

Fig1

DGN-WB joint heritability. Local h2 is estimated with SNPs within 1 Mb of each gene. Global h2 is estimated with SNPs that are eQTLs in the Framingham Heart Study on other chromosomes (FDR < 0.05).

otherfile<-my.dir %&% 'expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.all.Chr1-22_globalOtherChr.2015-03-18.txt'

fdrother<-read.table(otherfile,header=T) ##FHS eQTLs w/fdr<0.05 on non-gene chromosomes used to define global GRM

##Plot FDR based results
a<-ggplot(fdrother,aes(x=loc.jt.h2,y=glo.jt.h2)) + geom_point(cex=0.8) + geom_abline(intercept=1, slope=-1,color='red') + xlab(expression("local h"^2)) + ylab(expression("global h"^2)) + coord_cartesian(ylim=c(0,1),xlim=c(0,1)) + theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"))

##plot joint h2 estimates
local <- fdrother %>% select(loc.jt.h2,loc.jt.se) %>% arrange(loc.jt.h2) 
names(local) = c('h2','se')
data <- local %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE  TRUE 
##  4974  5989
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE  TRUE 
##  45.4  54.6
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.13
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))

my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05,  y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05,  y=0.66, hjust=0,gp=gpar(fontsize=14)))

b<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+coord_cartesian(ylim=c(0,1))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)

global <- fdrother %>% select(glo.jt.h2,glo.jt.se) %>% arrange(glo.jt.h2) 
names(global) = c('h2','se')
data <- global %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE  TRUE 
## 10505   458
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE  TRUE 
##  95.8   4.2
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.076
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))

my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05,  y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05,  y=0.66, hjust=0,gp=gpar(fontsize=14)))

c<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("global h"^2)) + xlab(expression("genes ordered by global h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)

multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)

tiff(filename=fig.dir %&% "Fig1.tiff",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "Fig1.png",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf 
##   2

Fig2

DGN-WB joint heritability with known trans-eQTLs. Local h2 is estimated with SNPs within 1 Mb of each gene. Known trans h2 is estimated with SNPs that are trans-eQTLs in the Framingham Heart Study for each gene (FDR < 0.05).

transfile<-my.dir %&% 'expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.all.Chr1-22_transForGene.2015-03-20.txt'

fdrtrans<-read.table(transfile,header=T) ##FHS trans-eQTLs for gene w/fdr<0.05 used to define Known trans GRM

##Plot FDR based results
a<-ggplot(fdrtrans,aes(x=loc.jt.h2,y=trans.jt.h2)) + geom_point(cex=0.8) + geom_abline(intercept=1, slope=-1,color='red') + xlab(expression("local h"^2)) + ylab(expression("known trans h"^2)) + coord_cartesian(ylim=c(0,1),xlim=c(0,1)) + theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"))
##plot joint h2 estimates
local <- fdrtrans %>% select(loc.jt.h2,loc.jt.se) %>% arrange(loc.jt.h2) 
names(local) = c('h2','se')
data <- local %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE  TRUE 
##   462   732
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE  TRUE 
##  38.7  61.3
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.133
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))

my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05,  y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05,  y=0.66, hjust=0,gp=gpar(fontsize=14)))

b<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax,color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+coord_cartesian(ylim=c(0,1),xlim=c(-30,1280))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)

global <- fdrtrans %>% select(trans.jt.h2,trans.jt.se) %>% arrange(trans.jt.h2) 
names(global) = c('h2','se')
data <- global %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE  TRUE 
##  1144    50
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE  TRUE 
##  95.8   4.2
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.021
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))

my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05,  y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05,  y=0.66, hjust=0,gp=gpar(fontsize=14)))

c<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax,color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("known trans h"^2)) + xlab(expression("genes ordered by known trans h"^2))+coord_cartesian(ylim=c(0,1),xlim=c(-30,1280))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)

multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)

tiff(filename=fig.dir %&% "Fig2.tiff",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "Fig2.png",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf 
##   2

Fig3

Polygenic v. sparse by elastic net.

data<-read.table(my.dir %&% 'working_DGN-WB_exp_10-foldCV_1-reps_elasticNet_eachAlphaR2_hapmap2snps_chr22_2015-01-21.txt',header=T,check.names=F)
ngenes<-dim(data)[1]
print("Elastic Net DGN-WB chr22 (" %&% ngenes %&% " genes)")
## [1] "Elastic Net DGN-WB chr22 (341 genes)"
data_long<-melt(data,by=gene)
## Using gene as id variables
a <- ggplot(data_long, aes(x = as.numeric(levels(variable))[variable] , y = value), group=gene) + geom_line(lwd=0.5,show_guide = FALSE,linetype=1) + aes(color = gene) + xlab(expression(paste("elastic net mixing parameter (",alpha, ")"))) + ylab(expression(paste("10-fold cross-validation R"^2))) + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))+ geom_point(show_guide = FALSE)
print(a)

b<-ggplot(data,aes(x=data[,12],y=data[,2]))+geom_point()+ ylab(expression(paste("R"^2," (",alpha," = 0)"))) + xlab(expression(paste("R"^2," (",alpha," = 0.5)"))) + geom_abline(intercept=0, slope=1,color='red') + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))

c<-ggplot(data,aes(x=data[,12],y=data[,22]))+geom_point()+ ylab(expression(paste("R"^2," (",alpha," = 1)"))) + xlab(expression(paste("R"^2," (",alpha," = 0.5)"))) + geom_abline(intercept=0, slope=1,color='red') + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))

data_pairs <- cbind(data[,3],data[,12],data[,21],data[,22])
colnames(data_pairs)<-c("alpha_0.05","alpha_0.5","alpha_0.95","alpha_1")
ggpairs(data_pairs)

multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)

tiff(filename=fig.dir %&% "Fig3.tiff",width=960,height=320)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "Fig3.png",width=960,height=320)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)
dev.off()
## pdf 
##   2

Fig4

GTEx tissue-wide heritability. Local h2 is estimated with SNPs within 1 Mb of each gene.

h2TW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=T)
seTW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_se.txt",header=T)

###make ordered point+CI h2 plots
gTW_h2 <- gather(h2TW,"Tissue.h2","TissueWide.h2",2:11)
gTW_se <- gather(seTW,"Tissue.se","TissueWide.se",2:11)
gTW_h2_se <- cbind(gTW_h2,gTW_se[,3])
colnames(gTW_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTW_h2_se)/length(unique(gTW_h2_se$Tissue))
gTW_h2_se <- gTW_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0) 
fig4<-ggplot(gTW_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),strip.text.x = element_text(size = 18),legend.text=element_text(size=14),legend.title=element_text(size=14))

###calc % nonzero for each tissue TISSUE-WIDE
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-gTW_h2_se %>% select(ensid,Tissue,nonzeroCI) %>% spread(Tissue,nonzeroCI)
for(i in 2:11){
  tis<-colnames(a)[i]
  cat("\n\n---",tis,"---\n")
  print(table(a[,i]))
  per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
  print(per)
  pvec <- c(pvec,per[2])
}
## 
## 
## --- CrossTissue ---
## 
## FALSE  TRUE 
## 14069  2938 
## 
## FALSE  TRUE 
##  82.7  17.3 
## 
## 
## --- AdiposeSubcutaneous ---
## 
## FALSE  TRUE 
## 15880  1127 
## 
## FALSE  TRUE 
## 93.40  6.63 
## 
## 
## --- ArteryTibial ---
## 
## FALSE  TRUE 
## 15836  1171 
## 
## FALSE  TRUE 
## 93.10  6.89 
## 
## 
## --- HeartLeftVentricle ---
## 
## FALSE  TRUE 
## 16258   749 
## 
## FALSE  TRUE 
##  95.6   4.4 
## 
## 
## --- Lung ---
## 
## FALSE  TRUE 
## 16078   929 
## 
## FALSE  TRUE 
## 94.50  5.46 
## 
## 
## --- MuscleSkeletal ---
## 
## FALSE  TRUE 
## 15944  1063 
## 
## FALSE  TRUE 
## 93.70  6.25 
## 
## 
## --- NerveTibial ---
## 
## FALSE  TRUE 
## 15541  1466 
## 
## FALSE  TRUE 
## 91.40  8.62 
## 
## 
## --- SkinSunExposedLowerleg ---
## 
## FALSE  TRUE 
## 15809  1198 
## 
## FALSE  TRUE 
## 93.00  7.04 
## 
## 
## --- Thyroid ---
## 
## FALSE  TRUE 
## 15680  1327 
## 
## FALSE  TRUE 
##  92.2   7.8 
## 
## 
## --- WholeBlood ---
## 
## FALSE  TRUE 
## 16062   945 
## 
## FALSE  TRUE 
## 94.40  5.56
###calc mean h2 for each tissue
a<-gTW_h2_se %>% select(ensid,Tissue,h2) %>% spread(Tissue,h2)
for(i in 2:11){
  tis<-colnames(a)[i]
  meanh2 <- signif(mean(a[,i]),3)
  cat(tis,"mean h2:",meanh2,"\n")
  h2vec <- c(h2vec,meanh2)
}
## CrossTissue mean h2: 0.057 
## AdiposeSubcutaneous mean h2: 0.0338 
## ArteryTibial mean h2: 0.0358 
## HeartLeftVentricle mean h2: 0.0383 
## Lung mean h2: 0.0325 
## MuscleSkeletal mean h2: 0.0279 
## NerveTibial mean h2: 0.044 
## SkinSunExposedLowerleg mean h2: 0.0351 
## Thyroid mean h2: 0.0392 
## WholeBlood mean h2: 0.0265
ann_text <- data.frame( h2 = rep(0.78,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
fig4.2<-fig4+geom_text(data=ann_text,aes(label=paste("mean_h^2 ==",mean_h2,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5) 
ann_text <- data.frame( h2 = rep(0.9,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
fig4<-fig4.2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)

fig4

tiff(filename=fig.dir %&% "Fig4.tiff",width=720,height=960)
fig4
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "Fig4.png",width=720,height=960)
fig4
dev.off()
## pdf 
##   2

Fig5

GTEx cross-tissue and tissue-wide h2 (A) and SE (B).

h2TW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=T)
seTW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_se.txt",header=T)

gh2_TW<-gather(h2TW,"CrossTissue","Tissue",3:11)
colnames(gh2_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5a<-ggplot(gh2_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Wide h'^2)) +  ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) + theme(plot.title = element_text(hjust = 0,face="bold")) 

gse_TW<-gather(seTW,"CrossTissue","Tissue",3:11)
colnames(gse_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5b<-ggplot(gse_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab('Cross-Tissue SE') + xlab('Tissue-Wide SE') +  ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+ theme(plot.title = element_text(hjust = 0,face="bold")) 
multiplot(fig5a,fig5b)

tiff(filename=fig.dir %&% "Fig5.tiff",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "Fig5.png",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## pdf 
##   2
chr22<-read.table(my.dir %&% "gtex-annot/gencode.v18.genes.patched_contigs.summary.protein.chr22")
gtexEN<-read.table(my.dir %&% "GTEx_tissue-wide_elasticNet_for_ggplot2_2015-05-15.txt",header=T)
chr22EN <- filter(gtexEN,ensid %in% chr22$V5)
ngenes <- length(unique(chr22EN$ensid))
ggplot(chr22EN, aes(x = alpha , y = R2), group=ensid) + facet_wrap(~tissue,ncol=3) + geom_point(show_guide = FALSE) + geom_line(lwd = 0.5, show_guide = FALSE) + aes(color = ensid) + xlab("alpha") + ylab("R2") + ggtitle("GTEx tissue-wide (" %&% ngenes %&% " chr22 genes)")

topR2<-filter(gtexEN,R2>0.5)
ngenes <- length(unique(topR2$ensid))
ggplot(topR2, aes(x = alpha , y = R2), group=ensid) + facet_wrap(~tissue,ncol=3) + geom_point(show_guide = FALSE) + geom_line(lwd = 0.5, show_guide = FALSE) + aes(color = ensid) + xlab("alpha") + ylab("R2") + ggtitle("GTEx tissue-wide (" %&% ngenes %&% " genes with R2>0.5)")

FigS1

GTEx tissue-specific heritability. Local h2 is estimated with SNPs within 1 Mb of each gene.

h2TS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=T)
seTS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_se.txt",header=T)

###make ordered point+CI h2 plots
gTS_h2 <- gather(h2TS,"Tissue.h2","TissueSpecific.h2",2:11)
gTS_se <- gather(seTS,"Tissue.se","TissueSpecific.se",2:11)
gTS_h2_se <- cbind(gTS_h2,gTS_se[,3])
colnames(gTS_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTS_h2_se)/length(unique(gTS_h2_se$Tissue))
gTS_h2_se <- gTS_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0) 
figS1<-ggplot(gTS_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),strip.text.x = element_text(size = 18),legend.text=element_text(size=14),legend.title=element_text(size=14))+ coord_cartesian(ylim=c(0,1))

###calc % nonzero for each tissue TISSUE-WIDE
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-gTS_h2_se %>% select(ensid,Tissue,nonzeroCI) %>% spread(Tissue,nonzeroCI)
for(i in 2:11){
  tis<-colnames(a)[i]
  cat("\n\n---",tis,"---\n")
  print(table(a[,i]))
  per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
  print(per)
  pvec <- c(pvec,per[2])
}
## 
## 
## --- CrossTissue ---
## 
## FALSE  TRUE 
## 14069  2938 
## 
## FALSE  TRUE 
##  82.7  17.3 
## 
## 
## --- AdiposeSubcutaneous ---
## 
## FALSE  TRUE 
## 16800   207 
## 
## FALSE  TRUE 
## 98.80  1.22 
## 
## 
## --- ArteryTibial ---
## 
## FALSE  TRUE 
## 16701   306 
## 
## FALSE  TRUE 
##  98.2   1.8 
## 
## 
## --- HeartLeftVentricle ---
## 
## FALSE  TRUE 
## 16709   298 
## 
## FALSE  TRUE 
## 98.20  1.75 
## 
## 
## --- Lung ---
## 
## FALSE  TRUE 
## 16784   223 
## 
## FALSE  TRUE 
## 98.70  1.31 
## 
## 
## --- MuscleSkeletal ---
## 
## FALSE  TRUE 
## 16573   434 
## 
## FALSE  TRUE 
## 97.40  2.55 
## 
## 
## --- NerveTibial ---
## 
## FALSE  TRUE 
## 16628   379 
## 
## FALSE  TRUE 
## 97.80  2.23 
## 
## 
## --- SkinSunExposedLowerleg ---
## 
## FALSE  TRUE 
## 16558   449 
## 
## FALSE  TRUE 
## 97.40  2.64 
## 
## 
## --- Thyroid ---
## 
## FALSE  TRUE 
## 16565   442 
## 
## FALSE  TRUE 
##  97.4   2.6 
## 
## 
## --- WholeBlood ---
## 
## FALSE  TRUE 
## 16517   490 
## 
## FALSE  TRUE 
## 97.10  2.88
###calc mean h2 for each tissue
a<-gTS_h2_se %>% select(ensid,Tissue,h2) %>% spread(Tissue,h2)
for(i in 2:11){
  tis<-colnames(a)[i]
  meanh2 <- signif(mean(a[,i]),3)
  cat(tis,"mean h2:",meanh2,"\n")
  h2vec <- c(h2vec,meanh2)
}
## CrossTissue mean h2: 0.057 
## AdiposeSubcutaneous mean h2: 0.0204 
## ArteryTibial mean h2: 0.0249 
## HeartLeftVentricle mean h2: 0.0323 
## Lung mean h2: 0.0223 
## MuscleSkeletal mean h2: 0.0214 
## NerveTibial mean h2: 0.0281 
## SkinSunExposedLowerleg mean h2: 0.0288 
## Thyroid mean h2: 0.0265 
## WholeBlood mean h2: 0.0261
ann_text <- data.frame( h2 = rep(0.78,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
figS1.2<-figS1+geom_text(data=ann_text,aes(label=paste("mean_h^2 ==",mean_h2,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5) 
ann_text <- data.frame( h2 = rep(0.9,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
figS1<-figS1.2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)

figS1

tiff(filename=fig.dir %&% "FigS1.tiff",width=720,height=960)
figS1
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "FigS1.png",width=720,height=960)
figS1
dev.off()
## pdf 
##   2

FigS2

GTEx cross-tissue and tissue-specific h2 (A) and SE (B).

h2TS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=T)
seTS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_se.txt",header=T)

gh2_TS<-gather(h2TS,"CrossTissue","Tissue",3:11)
colnames(gh2_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2a<-ggplot(gh2_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Specific h'^2)) +  ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) + theme(plot.title = element_text(hjust = 0,face="bold"))

gse_TS<-gather(seTS,"CrossTissue","Tissue",3:11)
colnames(gse_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2b<-ggplot(gse_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab('Cross-Tissue SE') + xlab('Tissue-Specific SE') +  ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+ theme(plot.title = element_text(hjust = 0,face="bold")) 
multiplot(figS2a,figS2b)

tiff(filename=fig.dir %&% "FigS2.tiff",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## pdf 
##   2
png(filename=fig.dir %&% "FigS2.png",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## pdf 
##   2